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When each site of a spatially extended excitable medium is independently driven by a Poisson 
stimulus with rate h, the interplay between creation and annihilation of excitable waves leads to 
an average activity F. It has recently been suggested that in the low-stimulus regime (h ~ 0) 
the response function F(h) of hypercubic deterministic systems behaves as a power law, F ~ h m . 
Moreover, the response exponent m has been predicted to depend only on the dimensionality d of 
the lattice, m = 1/(1 + d) [T. Ohta and T. Yoshimura, Physica D 205, 189 (2005)]. In order to 
test this prediction, we study the response function of excitable lattices modeled by either coupled 
Morris-Lecar equations or Greenberg-Hastings cellular automata. We show that the prediction is 
verified in our model systems for d = 1, 2, and 3, provided that a minimum set of conditions 
is satisfied. Under these conditions, the dynamic range — which measures the range of stimulus 
intensities that can be coded by the network activity — increases with the dimensionality d of the 
network. The power law scenario breaks down, however, if the system can exhibit self-sustained 
activity (spiral waves) . In this case, we recover a scenario that is common to probabilistic excitable 
media: as a function of the conductance coupling G among the excitable elements, the dynamic 
range is maximized precisely at the critical value G c above which self-sustained activity becomes 
stable. We discuss the implications of these results in the context of neural coding. 

PACS numbers: 87.19.L-, 87.10.-c, 87.19.1q, 87.18.Vf, 05.45.-a 



I. INTRODUCTION 

Sensory stimuli impinge continuously onto the periph- 
eral nervous system, where they are transduced into elec- 
trical activity of sensory neurons. Understanding how 
those and subsequent neurons encode and process the in- 
formation of the stimulus remains a formidable challenge 
for neuroscience since the pioneering work of Adrian [l|, 
and is the subject of ongoing research (see, e.g., Ref. [2[ 
for recent progress on olfaction). 

One of the most remarkable achievements of the ner- 
vous systems of multicellular organisms is their large dy- 
namic range, i.e., their ability to cope with stimulus in- 
tensities which vary by many orders of magnitude. Ex- 
perimental evidence in this direction is abundant, the 
simplest example being the century-old psychophysical 
laws: the psychological perception F of a given stimulus 
intensity h has been shown to be a power law for weak 
stimuli, F ~ h m . This behavior of the response curve 
F(h) is known as Stevens' law, and the response expo- 
nent m is called Stevens' exponent in the psychophysical 
literature [|[. Microscopic (i.e., neural) data also confirm 
this scenario: the activity of relay stages in sensory pro- 
cessing also increases as a power law of the stimulus inten- 
sities (e.g. glomeruli and mitral cells for olfaction [1,0, 
or ganglion cells of the retina [1, 0] ) • hi both cases (psy- 
chophysical and neural) , the response exponents are typ- 
ically less than 1, which indicates (as we will see below) 
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a large dynamic range of the response curves. 

That large dynamic ranges should be evolutionarily fa- 
vorable is generally agreed upon, owing to the fact that 
natural stimuli indeed span several decades of intensity. 
However, experimental results show that the dynamic 
range of the very first sensory neurons which perform 
the initial transduction is usually small, their firing rate 
varying essentially linearly with stimulus intensity (see, 
e.g., Ref.[|[ for the case of olfaction). Therefore, what 
remains to be explained is how those apparently conflict- 
ing results can be reconciled. In other words, how can 
large dynamic ranges be implemented by neurons? 

Two main mechanisms have long been proposed. The 
first one is adaptation, by which neurons manage to ad- 
just their range of operation according to the statistics 
of the ambient stimulus 0, [l(| HH, Ell HH • The second 
one is the intrinsic variation of firing thresholds among a 
population of sensory neurons, which would allow them 
to cover a wide range of stimuli (in spite of each of them 
having a small range) [l4[ ■ Both mechanisms can indeed 
contribute to an enhancement of dynamic range. How- 
ever, note that neither adaptation nor threshold variation 
requires interactions among neurons to work, insofar as 
adaptation has been understood as a dynamical process 
which neurons undergo individually and the firing thresh- 
old of a sensory neuron in principle does not depend on 
the activity of other sensory neurons. Therefore, if these 
were the only mechanisms responsible for enhancement 
of sensitivity and dynamic range, there should be no sig- 
nificant change in those properties if lateral connections 
among neurons were blocked. 

Experimental data, however, suggest otherwise. Deans 
et al. Q have measured the response function (firing rate 
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vs light intensity) of retinal ganglion cells of mice. For 
wild- type mice, they found a class of cells that responded 
with large dynamic range. When the same experiment 
was repeated with connexin36 knockout mice (i.e., mice 
that lack electrical synapses), they found that both sen- 
sitivity and dynamic range were significantly reduced. 
This suggests a third mechanism for dynamic range en- 
hancement, based on the interaction among neurons. 

This third mechanism is the subject of the present con- 
tribution. Previous work has revealed that, when ex- 
citable neurons are coupled (via chemical or electrical 
synapses), the response function of the resulting excitable 
medium indeed has much enhanced sensitivity and dy- 
namic range @, El El El El El El El, as compared 
to those of isolated neurons. The underlying mechanism 
relies on very general properties of excitable media: in- 
coming stimuli generate excitable waves which will dis- 
appear (due to the nonlinearity of their dynamics) upon 
collision with one another and/or with the system bound- 
aries. For weak stimuli, waves are rare and can propagate 
a long way before annihilation, therefore amplification is 
large (as compared with what would be observed for un- 
coupled neurons); for strong stimuli, waves contribute lit- 
tle to the overall network activity (since most neurons are 
being externally driven), therefore amplification is small. 
As a result, the medium as a whole has much larger sensi- 
tivity and enhanced dynamic range as compared to those 
of its building blocks @, El El El El El El El ■ 

The above reasoning has been tested and confirmed 
in a variety of models. In Refs. [H, El El El the cou- 
pling among excitable elements was probabilistic (say, via 
a transmission rate A). In such a scenario, low-stimulus 
amplification as described above occurs via stochastic ex- 
citable waves, whose (finite) lifetimes are essentially pro- 
portional to A (for small A). The dynamic range then ini- 
tially increases with increasing A, up to a critical value 
A = A c , where the system undergoes a nonequilibrium 
phase transition. Above A c self-sustained activity be- 
comes stable (i.e., small fluctuations can lead to non-zero 
density of active sites even in the absence of external 
stimuli). This hinders the coding of weak stimuli (just 
as a whisper cannot be heard in a sound system domi- 
nated by audio feedback), a problem that only worsens 
if the coupling is further increased. The dynamic range 
then decreases above A c and one concludes that the max- 
imum dynamic range is obtained precisely at the phase 
transition (l8j . 

Due to their probabilistic nature, the above cited sys- 
tems were cast in a framework of stochastic lattice mod- 
els, from which useful insights could be obtained by ap- 
plying mean field approximations and relying on well- 
known results of the statistical physics of nonequilib- 
rium phase transitions. For instance, the response ex- 
ponent m at criticality was shown to be a critical ex- 
ponent [H, Eli I20I. l2lll whose value has been known for 
over two decades [22l|. This should be contrasted with the 
models employed in Refs. flU El, El , where the coupling 
among excitable elements was deterministic. In these pa- 



pers, the models were such that no self-sustained activity 
was observed for vanishing stimulus rates. Besides, even 
if a transition to the self-sustained regime occurred, the 
standard results from statistical physics would not be 
easily applicable due to the deterministic nature of the 
excitable waves. 

In this context, our aim here is to fill two gaps: first, we 
verify the existence of power law responses in determin- 
istic excitable media without self-sustained activity; sec- 
ond, we probe the robustness of these power laws. To ac- 
complish the first goal, we have chosen to simulate hyper- 
cubic excitable media. This allowed us to test a theoret- 
ical prediction which has recently been proposed (based 
on scaling arguments) for the dependence of the response 
exponent m on the dimensionality d (23j . Moreover, it 
reveals important differences (regarding the dependence 
of m on d) with systems where coupling is probabilistic 
(as recently studied .21]). To accomplish the second goal, 
we employed the same model to show that, with a small 
change in its parameters, self-sustained activity can oc- 
cur, thus setting limits on the validity of the theoretical 
prediction. As it turns out, this last result puts the de- 
terministic and probabilistic cases in a similar state of 
affairs, where the dynamic range is maximized precisely 
at the transition to self-sustained activity. 

This paper is organized as follows. In Sec. [Ill the two 
models employed are described. The response functions 
in the absence and presence of self-sustained activity are 
analyzed in Secs. lIII Al and lTlTBl respectively. From these 
response functions we obtain the dynamic range, which 
is dealt with in Sec. lIVl Our conclusions are summarized 
in Sec. El 



II. MODELS 

In our simulations, we make use of a lattice in which 
each excitable site i is governed by the Morris-Lecar (ML) 
equations [H,E3 
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where the membrane capacitance per unit area is C rn = 
1/iF/cm 2 , membrane voltages Vi are measured in mV, 
current densities in //A/cm 2 , 0=1/3 ms _1 , and max- 
imal conductances for calcium, potassium, and passive 
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membrane leakage are respectively Gc a — 1 mS/cm 2 , 
G K = 2 mS/cm 2 , and G m = 0.5 mS/cm 2 . The corre- 
sponding reversal potentials are Ec a = 100 mV, Ek = 
—70 mV, and V res t = — 35 mV. Note that the gating vari- 
able for calcium mca is assumed to be always in equilib- 
rium, while w (which gates potassium currents) obeys a 
first-order dynamics [2 51 ] (both are dimcnsionlcss). All 
times are expressed in milliseconds. 

Even though the ML equations were developed origi- 
nally to describe the membrane potential of the barnacle 
muscle fiber, our aim here is not to model any specific bio- 
logical tissue in particular, but rather to shed light on the 
influence of the network topology on its response prop- 
erties, particularly the dynamic range. Here we study 
hypcrcubic lattices with dimensionality d, restricting our- 
selves to the simplest case of electrical coupling, for which 
the synaptic currents are given by Ohm's law, 
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where j runs over the first neighbors of i. The conduc- 
tance Gij between sites i and j could account for gap 
junctions (e.g. as observed in axoaxonal contacts in the 
hippocampus [26|, [23] or dendrodendritic contacts of mi- 
tral cells in the olfactory glomeruli [28|) or ephaptic in- 
teractions (as modeled by Bokil et al. to occur in the 
olfactory nerve [29l|). 

The external current lf lm (t) accounts for the stimuli 
arriving in the network, which we model as a Poisson pro- 
cess. Each neuron independently receives current pulses 
at constant rate h (measured in ms" 1 ). Each pulse has 
duration D and intensity Iq (so that for h > D" 1 the 
regime of a continuous external current is approached). 

To test the robustness of the results and to allow for 
larger system sizes, we also simulate lattices in which each 
excitable element is modeled by the n-state deterministic 
Greenberg-Hastings cellular automaton (GHCA) [3(|. In 
this case, each site i at discrete time t can be in states 
Xi(t) G {0, 1, 2, . . . , n — 1}, where x = and 1 represent a 
quiescent (polarized) and spiking (depolarized) neuron, 
respectively, whereas for 2 < x < n — 1 the site is re- 
fractory The dynamical rules are cyclical: if Xi(t) > 1, 
then Xi(t + 1) = [xi(i) + 1] mod n, i.e., after a spike the 
model neuron deterministically undergoes n — 2 refrac- 
tory steps before returning to the x = quiescent state. 
If Xi(t) = 0, then Xi(t + 1) = 1 if at least one of its 
2d nearest neighbors is spiking at time t or if an exter- 
nal stimulus arrives at site i [xi(t + 1) = otherwise]. 
The Poissonian external stimulus occurs independently 
at each site with probability P = 1 — exp(— hd), where 
<5 = 1 ms is the time step adopted in this case. 

For both models, i — 1, . . . , N, where N = L d is the 
total number of excitable elements in a network of linear 
size L. 



III. RESPONSE OF HYPERCUBIC EXCITABLE 
MEDIA 

Let F be the mean firing rate, defined as the total 
number of spikes in an interval T max , divided by the 
number N of neurons and by T max . To avoid under- 
sampling in the low-stimulus regime, we have chosen 
T m .ax = max{n/(WV), 100 ms}, where n is the approx- 
imate mean number of attempts to initiate an excitable 
wave (we have typically employed n = 25). We define the 
response function (or transfer function) of the network to 
the external stimulus as F(h). In the following, we make 
use of a uniform coupling Gij = G and study how the 
response function changes with G. 



A. Power laws 

Figure Q] shows the results for a one-dimensional ML 
lattice with D = 0.3 ms and Iq = 15 /iA/cm 2 . As G in- 
creases, three regimes are observed in the response of the 
network. For weak coupling [left panel in Fig. trian- 
gles in Fig. [TJd] , synaptic currents from spiking neighbors 
are not strong enough to generate spikes, so each stimu- 
lus event generates one spike, and the response function 
increases linearly (up to saturation at F max , which is 
essentially the inverse of the refractory period). Above 
a certain value G[ ~ 0.14 mS/cm 2 , however, the con- 
ductance is strong enough to allow the propagation of 
excitable waves. In this regime [middle panel in Fig. [lji, 
squares in Fig. [TJd], which is observed up to a second 
transition at G'[ ~ 0.24 mS/cm 2 , excitable waves are 
created by external stimuli and annihilated by one an- 
other and by the boundaries (open boundary conditions 
have been employed throughout this paper). Above G'{, 
current leakage to neighbors is so large that it typically 
prevents neurons from spiking upon the incidence of a 
single stimulus pulse. What we observe [right panel of 
Fig. DJt] is that a neuron will fire only if it is at the 
boundary (in which case it has fewer neighbors and con- 
sequently less leakage) or if two stimulus pulses happen 
to arrive nearly consecutively (in a mimicry of temporal 
summation). Note that in the three panels in Fig. [T^, 
the seed of the pseudo-random-number generator is the 
same, so the spikes in the left panel coincide with stim- 
ulus pulses. In the right panel of Fig. [TJt, however, only 
the stimulus pulses that happened to fall right at the 
borders generated waves (all other visible perturbations 
are subthreshold, not spikes). In this regime, inevitably 
poor statistics ensues, except for large stimulus rates, as 
reflected in the G = 0.3 (circles) curve in Fig. [TJd. Note, 
however, that if a spike is finally produced, propagation 
of an excitable wave does occur, which explains why the 
response in this case is larger than for G < G[. 

The response curves in Fig. [T]d clearly show power laws 
F ~ h m in the low-stimulus regime. For G < G' 1: the re- 
sponse is linear (m = 1) and can be easily explained: for 
each stimulus pulse, a small number of spikes is gener- 
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FIG. 1: (Color online) (a) Membrane potentials of 100 ML 
neurons versus time with h = 10~ 2 ms" 1 for G = 0.1 (left 
panel), 0.2, (middle panel) and 0.3 mS/cm 2 (right panel). 
The seed of the pseudo-random number-generator is the same 
for the three values of G. (b) Response function for a one- 
dimensional lattice of L = 1000 ML excitable elements. Sym- 
bols (bars) represent averages (standard deviations) over 10 
runs. Solid lines are power laws discussed in the text. 



FIG. 2: (Color online) (a) Snapshots of networks with 50 x 50 
ML neurons, with depolarized (spiking) membrane poten- 
tials coded as white. From left to right, G — 0.1, 0.5, and 
0.9 mS/cm 2 , with h = 2 x 10~ 3 ms" 1 and t = 3.0 ms (3.5 ms) 
for top (bottom) row. The seed of the pseudo-random-number 
generator is the same for the three values of G. (b) Response 
function for a network of 200 x 200 ML excitable elements. 
Symbols (bars) represent averages (standard deviations) over 
ten runs. Solid lines are power laws discussed in the text. 



ated (typically one) and excitable waves do not interact. 
For G[ < G < G'l, however, excitable waves are created 
in randomly located points and annihilate upon encoun- 
tering one another. To understand how this nonlinear 
interaction leads to a power law in the dependence of 
F on h, Ohta and Yoshimura have recently proposed an 
elegant scaling reasoning (23|. In the scaling regime, F 
should depend on a dimcnsionlcss variable A. Since h 
is small, the characteristic times for wave creation and 
wave annihilation are much smaller than the time of free 
propagation. Therefore, the only relevant parameters are 
the width I of an excitation, the wave speed c, and the 
rate h. Recalling that h is measured in events per unit 
time per site (thus having dimension of we ob- 

tain A = /ic _1 ^ 1+d . If we now assume a scaling relation 
F ~ A m , the exponent m can be obtained by noting 
that in the low-stimulus regime waves are sparsely dis- 
tributed and the dependence of F on I must be linear; 



hence m = 1/(1 + d) [H]. 

As shown in Fig. [TJd, this prediction is confirmed in 
our one-dimensional ML simulations in the parameter re- 
gion where excitable waves propagate ballistically. Par- 
ticularly for d — 1, the scaling relation F ~ h 1 ! 2 had 
already been conclusively confirmed for the GHCA (in 
both simulations @, [H| and analytical calculations Q) 
and coupled map lattices [Tfij . However, for more realis- 
tic models, it was only approximatel y v erified for a chain 
of Hodgkin-Huxley model neurons [Toll and a reaction- 
diffusion partial differential equation [23| , with exponents 
around m ~ 0.4. In Fig.[T|D we fill this gap with an agree- 
ment over more than two decades. 

In two dimensions, simulations have been carried out 
with stronger stimulus pulses (D = 0.45 ms and Iq = 
150 /zA/cm 2 ) to prevent excessive leakage owing to the 
larger number of neighbors. The same scenario has been 



5 



observed. For small G, each stimulus pulse generates at 
most an evanescent wave with a radius of a few neighbors 
[left panel of Fig. Ufa,]. For G > G' 2 ~ 0.225 mS/cm 2 , how- 
ever, generated waves can propagate ballistically with 
their radii increasing indefinitely. As shown in the middle 
panel of Fig. [2ji, annihilation in this case is more compli- 
cated than for d — 1, for now colliding waves may have 
different radii and their surfaces merge to form irregular- 
shaped excitations [H, [ij], H3 • This regime breaks down 
for G = G'2 — 0.725 mS/cm 2 , above which current leak- 
age is again too strong and spikes are generated with 
at least two nearly consecutive stimulus pulses or at the 
boundaries. Note that several waves that appear in the 
middle panel of Fig. [2^ are absent in the right panel, as 
exemplified by the arrows (as in Fig. [IJ the seed is the 
same for the three panels). Several waves in the right 
panel of Fig. [2^ have been created at the borders (and 
propagate faster than those of the middle panel because 
G is larger). 

As for the response functions, Fig. [2)3 shows that Ohta 
and Yoshimura's exponent m = 1/3 for d = 2 agrees 
(for two decades) with simulations for G = 0.5 mS/cm 2 
(and this holds true in the whole interval G" 2 < G < 
G'^)- Interestingly, another exponent (not predicted orig- 
inally [23|) arises for G > G'{\ in this regime, waves 
typically require two nearly consecutive stimulus pulses 
to be created, and for weak stimuli this occurs approx- 
imately at a rate h! ~ h 2 . But once they are created, 
Ohta and Yoshimura's reasoning is still valid, now with 
the dimensionless variable rewritten as A = h'c~ 1 l 1+d . 
We therefore obtain the exponent m = 2/(1 + d), which 
is reasonably confirmed for G = 0.9 mS/cm 2 (circles) 
in Fig. Looking back to the analogous situation for 
the one-dimensional case, the circles in Fig. [1}d are com- 
patible with an exponent m = 1 (the extremely poor 
statistics notwithstanding). Whether further increasing 
G leads to other transitions (inducing the necessity of, 
say, k > 2 nearly consecutive pulses to generate a wave) 
and new exponents [presumably m = k/(l + d)] is a ques- 
tion beyond the scope of this work, but perhaps worth 
pursuing. It is important to point out, however, that 
these transitions may have limited biological applicabil- 
ity: chemical synapses (which do not suffer from leakage) 
are not included in this model, yet abound in the nervous 
system. 

In order to test Ohta and Yoshimura's prediction in 
three dimensions, we have performed simulations of the 
GHCA model. With the rules defined in Sec. HH an in- 
coming stimulus pulse generates an excitable wave which 
propagates ballistically until annihilation with another 
wave or with the system borders [l?], Hfil, H3| , precisely 
as observed in the intermediate region G' < G < G" 
for the ML equations. The motivation for switching to 
a simpler model is that it allowed us to simulate much 
larger networks than would be feasible with the ML equa- 
tions. As shown in the response functions of Fig. [31 
finite-size effects arc strong. However, for a network of 
N = 160 3 automata (a system size beyond our computa- 



tional resources for the ML equations), it is already pos- 
sible to verify the power law F ~ /1 1 / 4 for more than two 
decades. Incidentally, we note that the response function 
of two-dimensional GHCA networks has been studied in 
Ref. [l7|, but the power law has been missed. The inset 
of Fig. [3] confirms the predicted exponent. 



1 



cU3 


« h 1,4 _ 




/L=160 










d = 2 






1=40 / / 
/l=20 f 


10- 1 

L=80 

10- 2 








r / L=1 ° 


10~ 7 10~ 5 10" 3 10~ 1 


10 1 



3 LZ C I 1 i 1 1 i I 

10~ 7 10 6 10~ 5 10~ 4 10~ 3 10~ 2 10~ 1 10° 10 1 
h (ms~ 1 ) 



FIG. 3: Response curves of the GHCA model in d = 3 for in- 
creasing system size (n — 3). Inset: GHCA response function 
for d = 2 (L = 2744, averages over five runs). Dashed lines 
shows the response exponent m = 1/(1 + d) for both cases. 



B. Spiral waves 

What we have described so far suggests that the re- 
sponse exponent is indeed m = 1/(1 + d) whenever the 
following two conditions are satisfied: (A) every quies- 
cent neuron (i.e., not only those at the borders) spikes 
upon the incidence of a single stimulus pulse and (B) 
this spike creates a deterministic excitable wave which 
will be annihilated at the borders or upon encountering 
other wave(s). In the examples shown in figures [T] and [2l 
these two conditions are simultaneously satisfied only for 
G' < G < G" . For G < G', condition A is satisfied, but 
B is not; for G > G" , condition B is satisfied, but A is 
not. 

The above scenario, however, is not general. As has 
been known for many decades, excitable media can ex- 
hibit self-sustained activity in the form of spiral or scroll 
waves, a topic that has received much attention due to its 
relevance in different scientific branches such as cardiol- 
ogy [Hl,[12|i cytology (HI, physics [HI, chemistry [35ll36j|. 
and neuroscience [37| . among others. In Fig.[2]the param- 
eters of the ML system were such that spiral waves did 
not occur. With a slight deviation in parameter space, 
however, spiral waves may appear, even in a homoge- 
neous lattice. As shown in the right panel of Fig. [4]i, this 
is the case for <fi = 0.4 ms" 1 and G = 0.35 mS/cm 2 (all 
other parameters remaining the same), for instance. In 
this system, spiral waves emerge [sec arrow in Fig. be- 
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FIG. 4: (Color online) Spiral waves with <f> — 0.4 rns -1 . 
(a) Snapshots of networks with 50 x 50 ML neurons, with 
depolarized (spiking) membrane potentials coded as white. 
From left to right, G = 0.15, 0.25, and 0.35 mS/cm 2 , with 
h = 4 x 10 -3 rns -1 and t = 20 ms (21 ms) for top (bottom) 
row. The seed of the pseudo-random-number generator is the 
same for the three values of G. (b) Response function for a 
network of 200 x 200 ML excitable elements. Symbols repre- 
sent averages over five runs (standard deviations are smaller 
than symbol size). Firing rates were measured over 150 ms 
after a 50 ms transient. Horizontal and vertical lines illustrate 
the relevant quantities for calculating the dynamic range A 
(see Eq. (7) and text for details). 



cause of the local inhomogcneities created by the stochas- 
tic input [lH and, once established, they typically re- 
sist being destroyed by the same stochastic input (even 
though their shape is continuously perturbed by the Pois- 
son pulses). 

With this new phenomenon at play, how does the sce- 
nario evolve as the coupling G changes? For low G (say, 
G < G' 2 ), the overall behavior of the system is the same 
as that of Fig. i.e., a stimulus-induced spike at one 
site does not propagate too far [compare the left panels 
of Figs. ®l and 2k]. Correspondingly, the response func- 
tion is linear. If G is increased, a transition occurs which 
allows the wave radii to increase indefinitely [Fig. [2k. and 



Fig. [4k, middle panels]. However, contrary to what was 
previously observed, this dynamic regime is no longer 
valid in a broad range of G values. As G is further in- 
creased, spiral waves quickly emerge. As for the second 
transition previously observed at G' 2 ', it now essentially 
loses meaning, for as soon as the waves are created — no 
matter whether by one or two incoming stimuli, or at the 
boundaries — the conditions are set for the spiral waves 
to dominate the network. 

Regarding Ohta and Yoshimura's conjecture in this 
scenario, the response function near the transition to 
self-sustained activity suffers from strong statistical fluc- 
tuations, as expected [see solid squares in the inset of 
Fig. [SJi] . It seems compatible with a power law with ex- 
ponent m = 1/3, but for less than a decade only [note 
that even the self-sustained activity suffers from finite- 
size effects for low enough stimulus rate — see pentagons 
in the inset of Fig. [SJi] . It is at present unclear whether 
larger system sizes or longer stimulus times would con- 
firm the power law at the transition. 

The drastic consequences of this self-sustained activity 
for the response curve are shown in Fig. UJd: the weak- 
stimulus response no longer decreases as a power law for 
decreasing h, but reaches instead a value Fq which corre- 
sponds to the average firing rate when the lattice is dom- 
inated by spiral waves. To obtain a reasonable estimate 
of Fq, we simulated the following protocol: 150 x 150 net- 
works were stimulated during a period T st i m = 100 ms 
with a constant rate h = 4 x 10 -3 ms -1 . The stim- 
ulus was then switched off (h = 0) and the mean ac- 
tivity of the network Fq was measured after a transient 
Ttrans = 900 ms. Figure [St shows how Fq depends on 
G. A transition is clearly seen near G ~ 0.275 mS/cm 2 , 
above (below) which Fq > (Fq = 0). In the inset of Fig- 
ure [5]; we also show the probability p that spiral waves 
survive after the transient, which was estimated by di- 
viding the number of runs in which spiral waves survived 
by the total number of runs. The sharpness of the p(G) 
curve also suggests a transition to a regime where self- 
sustained activity is stable. 

This second scenario appears to be more general than 
the one described in Sec. IIII Al We have simulated 
networks in which each element was modeled by the 
Hodgkin-Huxley equations [39[ with standard parame- 
ters [40[ and have obtained spiral waves. Moreover, one of 
the most studied causes of spiral wave creation is disorder 
and noise in the excitable dynamics [U [42|, HH , which 
are absent from the present study. We have nonethe- 
less tested some ML networks where (f> was distributed 
around 1/3 ms^ 1 with some variance, and have again 
obtained spiral waves. It is important to remark that, 
for the purposes of the present study, it is not enough 
that an excitable medium be able to sustain spiral waves 
in the absence of stimulus, say, for a given initial con- 
dition. The question is whether the Poisson stimulus is 
able to create spiral waves and, at the same time, allow 
them to survive. Consider, for instance, the limit of very 
weak stimuli (h — > 0). In this regime spiral waves hardly 
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emerge (even for <f> = 0.4 ms" 1 ) because fluctuations are 
not sufficiently strong [see, e.g., the open pentagons in the 
inset of Fig. [5ji] - At the other extreme, a large value of 
h can easily provide the necessary fluctuations, but then 
the created spiral waves will be statistically overshad- 
owed by the very stimuli that generated them. Overall, 
the probability of self-sustained activity coexisting with 
the Poisson stimulus depends not only on the model pa- 
rameters (in this case, <j> or G) but also on the system size 
(N), stimulus rate (h), and duration (T max ). A more de- 
tailed study of this dependence would be welcome. 



IV. DYNAMIC RANGE 

We can now return to the quantity that originally mo- 
tivated this study. The dynamic range A of a response 
curve F(h) is formally defined as [44| 



10 log 



id 



^0.9 

h .i 



(7) 



where ft-o.i (^0.9) is the stimulus intensity such that the 
difference F — F is 10 % (90 %) percent of the response 
interval F max — Fq. As depicted in Fig. [4b>, A mea- 
sures (in decibels) the range of stimulus intensities that 
can be "appropriately" coded by the mean firing rate of 
the system, discarding intensities whose corresponding 
responses are too close either to saturation (h > ho.g) or 
to baseline (h < /10.1). This measure of appropriateness 
is evidently arbitrary, but standard in the biological liter- 
ature and very useful, since it is a dimensionless quantity 
that allows direct comparison with experimental results. 

Figure O shows the behavior of the dynamic range (es- 
timated numerically from the response curves) as a func- 
tion of the coupling conductance. For d = 1 [Fig. [3^], A 
changes very little for G < G^, staying in the range of 
16 dB (which is comparable to experimental values of iso- 
lated olfactory sensory neurons [8[ and retinal ganglion 
cells of connexin36 knockout mice [1, 0] ) ■ The transition 
near G\ seems abrupt, after which the dynamic range be- 
comes substantially larger: the system attains ~ 31 dB, 
an enhancement of about 100 % which had also been 
previously obtained with a cellular automaton model Q • 
This enhancement is clearly due to a change in the re- 
sponse exponent m, which greatly amplifies weak stimuli 
[recall the squares in Fig. [lb] . For G > G'{ the dynamic 
range is reduced, once more because of the change in the 
weak-stimulus sensitivity [recall the circles in Fig. [Us]. It 
is important to remark that the poor statistics in Fig.[Tb> 
do not compromise the accuracy of the measured dy- 
namic range, since the strong fluctuations occur below 
the sensitivity threshold /lo.i- 

For <fi = 1/3 ms" 1 , the results in d = 2 are similar [see 
Fig. [5b] . As G' 2 is approached from below, the transition 
is somewhat smoother than for d = 1. More importantly, 
since the response exponent for G' 2 < G < G 2 (m =1/3) 
is smaller than for the corresponding regime in d = 1 
(m = 1/2), the weak-stimulus amplification for d = 2 



is larger and so is the dynamic range, which reaches ~ 
38 dB. The same trend in the dependence of A on d 
is observed in the GHCA model: with a fixed system 
size N = 14 6 , by varying the dimensionality, we obtain 
A = 31, 43 and 54 dB for d = 1, 2, and 3, respectively. 
Note that these values are comparable to those obtained 
in the ML model (the differences in d = 2 being explained 
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by finite-size effects, as extensively discussed in Rcf. 

This picture changes qualitatively when spiral waves 
come into play (0 = 0.4 ms , rightmost column of 
Fig. [5]). For G < G' 2 the dynamic range increases mono- 
tonically with G, reaching a maximum near G' 2 . Increas- 
ing G further, however, leads to the onset of spiral waves, 
and the nonzero baseline activity Fq prevents the appro- 
priate coding of weak stimuli. This is clearly seen in 
Fig. Eb> (pentagons): an observer would have much dif- 
ficulty in distinguishing the responses of any two points 
below h = 10~ 3 ms" 1 , which leads to a drastic decrease 
in dynamic range. Moreover this problem becomes more 
and more severe as G is further increased: since Fq in- 
creases with G for G > G' 2 [see Fig. [5b], the dynamic 
range decreases with increasing G. Therefore, if a deter- 
ministic excitable medium supports spiral waves in some 
parameter region, its dynamic range will be maximum 
precisely at the transition where they become stable. 



V. CONCLUDING REMARKS 

We have simulated hypercubic networks of excitable 
elements modeled by the Morris-Lccar equations and 
Greenbcrg-Hasting cellular automata. We have studied 
how the collective response F of the network to a Pois- 
son stimulus with rate h changes with the coupling G and 
the dimensionality d. Two scenarios have been observed. 
In the first one, a broad range of G values exists such 
that excitable waves are created and thereafter propa- 
gate ballistically, being annihilated upon encountering 
one another or the system boundaries. In this regime, 
the response function F(h; d) is shown to be a power law 
F ~ h m . Furthermore, we have confirmed that, if waves 
are created upon the incidence of a single stimulus pulse, 
the response exponent agrees with the theoretical predic- 
tion of Ohta and Yoshimura, m = 1/(1 + d) [23| . We have 
argued that in a regime where wave creation requires the 
incidence of two nearly consecutive stimuli, an exponent 
m = 2/(1 + d) should be expected and is confirmed by 
our ML simulations in d = 2 (also for a broad range of 
G values). 

If a system is such that the exponent m = 1/(1 + d) 
holds, the dynamic range increases with the dimension- 
ality d (as confirmed here for d = 1 and 2 in the ML 
model and d = 1, 2, and 3 for the GHCA model). This 
is in stark contrast with probabilistic excitable systems, 
where the maximum dynamic range attained at a given 
dimension d is a decreasing function of d. This hap- 
pens because in that case m corresponds to the critical 
exponent S^ 1 (apparently belonging to the directed per- 
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FIG. 5: (Color online) Left (right) column: absence (emergence) of spiral waves. Dynamic range versus coupling conductance 
for (f> = 1/3 ms" 1 [(a) and (b)] and <j> = 0.4 ms -1 (d). Triangles denote G < G' and squares denote G' < G < G" . Circles 
denote G > G" in the absence of self-sustained activity (b), whereas pentagons denote the spiral wave regime (d). In (c), 
the self-sustained activity Fo in the absence of stimulus is plotted against G for fixed <j> = 0.4 ms -1 . Inset of (c): estimated 
probability p of spiral wave survival versus G (see text for details). Inset of (d): response functions with <j> — 0.4 ms -1 for 
G = 0.25 (solid squares) and 0.275 mS/cm 2 (open pentagons). System sizes are L = 1000 in (a); N = 200 2 in (b) and (d); and 
N = 150 2 in (c). 



eolation universality class [HI), and S h 1 increases with 
d. 

In this context, one should not be misled by the appar- 
ent paradox posed by the assumption that a determin- 
istic system is "just" a particular case of a probabilistic 
one. Consider, for instance, a probabilistic version of 
the d-dimcnsional GHCA in which a stimulus would be 
transmitted to its quiescent neighbors with probability 
q: the function A(q) has qualitatively the same shape as 
that of Fig. [5ji and the maximum value of A attained at 
given d is a decreasing function of d i2lj. Why then for 
q = 1 do we have an increasing A(d)? Remember that 
the main condition for the exponent m = 1/(1 + d) to 
hold is the absence of self-sustained activity. In a prob- 
abilistic system, this requires not only that q is precisely 
1, but also that the initial conditions are appropriately 
set [13, H3 ■ For q infinitesimally smaller than 1 or q = 1 
with random initial conditions, self-sustained activity en- 
sues in the probabilistic GHCA. Therefore, in this partic- 
ular model the result m = l/(l+d) is obtained only under 
very artificial circumstances, at the edge of the param- 
eter space and only for restricted initial conditions. In 
contrast, for the deterministic ML lattices studied here, 
the exponent holds in a broad region of the parameter 



space for any initial condition. 

A substantially different scenario has been obtained 
with a change in a single parameter of the ML model, for 
which stable spiral waves were observed when the cou- 
pling was increased above a certain critical value (lead- 
ing to a breakdown of Ohta and Yoshimura's prediction). 
Given the ubiquity of spiral waves in studies of excitable 
media, this scenario is likely to be more general than the 
one previously described. In this case, a unifying picture 
emerges for both deterministic and probabilistic excitable 
media: the dynamic range in both cases is maximized at 
the critical value of coupling above which self-sustained 
activity becomes stable. 

Put into a broader context, our results reinforce the 
idea that optimal information processing near critical- 
ity, a topic which has received much attention in recent 
decades [45| , could have a bearing on the brain sciences. 
In fact, experimental results that are consistent with the 
hypothesis of neurons collectively operating near a criti- 
cal regime have recently appeared 0, [i?], HH, HU , joined 
by theoretical efforts aimed at understanding the compu- 
tations themselves [TH, [5(3, HH, [H| as wei l as t ne homc- 
ostatic mechanisms that could maintain the system at 
criticality [H, [54|. These issues still pose remarkable 
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challenges for the years to come, which opens the possi- 
bility of new lines of research connecting physicists with 
systems biology in general, and neuroscicncc in particu- 
lar. 
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